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ABSTR.^CT 



The vertical discretization in a linearized baroclinic prediction model was 
analyzed by comparing various finite element and finite difference solutions following 
Jordan (1985). Modifications were made on Jordan's (1985) Galerkin finite element 
approximation for two staggered grids. Comparisons were made with the unmodified 
models (Jordan, 1985) and with finite difference approximations for the same two 
staggered grids. The models were run with four experiments. Most of the oscillations 
that occurred in the temperature profiles near the surface of the unmodified Galerkin 
finite element approximations disappeared following the modifications. 
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I. INTRODUCTION 



Most numerical weather prediction models use finite differences to accomplish 
the vertical discretization even though they use fmite difference, finite element, or 
spectral, horizontal discretizations. The only exceptions are the Canadian regional and 
hemispheric models (Staniforth and Daley, 1977 and 1979), which use finite elements in 
the vertical. The successful numerical prediction of synoptic evolutions requires a 
proper representation of the vertical variation of the predictive fields. Since smaller 
scale features such as fronts (Hoskins and Bretherton, 1972 and Williams, 1967) and 
the large scale planetar>' waves (Gall, 1976) are forced by energetic synoptic-scale 
features, it follows that all predictive scales of motion may be sensitive to the vertical 
discretization used in the numerical models. 

Most of the finite difference vertical discretizations use a staggered arrangement 
of variables. It has been demonstrated that staggering of variables in the horizontal 
(Winninghoff, 1968; Arakawa and Lamb, 1977 and Schoenstadt, 1980) improves 
geostrophic adjustment and the response to small scale forcing. Most quasi-geostrophic 
models (Charney and Phillips, 1953) use vertical staggering where the vertical motion 
and the temperature are carried between the levels which carry horizontal velocity and 
pressure. This arrangement will be referred to as grid B. Lorenz (1960) introduced a 
different grid for the balance equations which was designed to conserve energy. This 
arrangement places only the vertical velocity between the levels which carry the other 
variables (horizontal velocity, pressure and temperature) and will be referred to as grid 
A. Tokioka (1978) analyzed a number of vertical grids with linearized equations and 
found that grid A has a computational mode in the temperature field. Arakawa (1984) 
compared baroclinic instability for grids A and B in the linearized quasi-geostrophic 
equations. He found a false short wave instability for grid A which did not occur with 
grid B. This problem is related to the computational mode in the temperature field. 
Another difficulty with grid B is that the matrix which must be inverted to find the 
temperature from the pressure is singular. This is especially important for initialization. 
Many operational primitive equation models use grid A for energy conservation. 

The use of finite elements for the vertical discretization can be expected to give a 
more accurate representation of vertical variations. The finite element method is a 
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special case of the Galerkin procedure which represents the dependent variables with a 
weighted sum of basis functions that have a prescribed spatial structure. The finite 
element method employs basis functions which are zero except in a limited region 
w'here they are low’-order polynomials. This method w'as developed in engineering 
statics (see e.g., Zienkiewhcz, 1977) and it has been more recently applied to fluid 
dynamics and hydrology (see Gray and Finder, 1976). The finite element method has 
been successfully applied to meteorological prediction with the shallow water equations 
by Cullen (1973), Hinsman (1975) and Staniforth and Mitchell {1977, 1978). Cullen 
(1973), Neta et al (1986), and Neta and Williams (1986) demonstrated that finite 
element formulations with piecewise linear basis functions are more accurate than 
second order finite differences. 

Jordan (1985) compared six linear, baroclinic, vorticity-divergence equation 
models using three grid schemes, grid A, grid B and an unstaggered grid. A 
perturbation w’as found in the temperature fields of the Rossby wave experiment and 
the mountain topography experiment in the finite element models using grids .A. and B 
(Fig. 1.1). 

The purpose of this study is to see if the perturbations noted in the Jordan study 
could be fixed. A baroclinic instability experiment is made with a comparison of finite 
element and finite difference models for grids A and B. The finite element models for 
grids A and B are modified at the boundaries and the finite difference models for grids 
A and B are left as they were except for the necessary modifications to run the 
baroclinic instability experiment. The results of the modified models are compared 
with the results of the unmodified models from Jordan (1985) for the Rossby wave 
experiment, the mountain topography experiment and the diabatic heating experiment. 
The experiments are described in Chapter III. 
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Fig. 1.1 Two vertical grids. 



12 



II. MODEL DESCRIPTIONS 



A. MODEL FEATURES 

Jordan (1985) developed six numerical models with several features to make easy 
modifications for a wide range of experiments. A menu is added to each of the models 
to make the transition between four experiments simple. The user is able to prescribe 
heating, mountain topography, velocity perturbation, or baroclinic experiments and the 
model will make the prescribed changes in the variables governing these cases. Another 
menu is added so that the user can prescribe the amount of printout desired from each 
run. The models are written in modular structure using FORTRAN 'll. There is 
parallel construction between models. The subroutines used in one model are very 
similar to those used in the other models. The models run quickly on an IBM-3033 
mainframe; for example, a 96-hour forecast for a 12 layer finite element model uses less 
than five seconds of computer processing time. 

B. GOVERNING EQUATIONS 

Each model approximates the same set of governing equations. The vorticity 
equation (2.1), the divergence equation (2.2), the surface geopotential equation (2.3), 
and the first law of thermodynamics (2.4) are the prognostic equations for the forecast 
variables vorticity, divergence, surface geopotential and potential temperature. The 
surface geopotential equation is the lower boundary condition on the vertical velocity. 
The vertical coordinate Z = — ln(p/pQ) is used, but the non-Boussinesq terms 

involving e*^ are replaced by one. The prognostic equations in the coordinates x, y, Z, 
and t are: 



dC ^ 5w dw 5w 5u 

— + (^ + f)D + pv + = 0, 

dt ^ ^ dx. dZ dydZ 



( 2 . 1 ) 



dD -f- 

dt ^x dx dy dy dy dx dx dZ 
dw dv ^ 

+— — + pu - fi; + 0^ 

dy dZ 



( 2 . 2 ) 
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(2.3) 



d(p5 

— ^ = MTS , 
dt 



dT 



(2.4) 



In these equations: 

q is the vertical component of vorticity, ~ dvldx — 5u, ^y, 

D is the horizontal divergence, D = ou/^x + dv/dy, 

(p is the geopotential, (p = gZ, 

(pj is the surface geopotential, 

T is the potential temperature, 
u is the x-component of velocity, 

V is the y-component of velocity, 
w is the vertical velocity, 

Q is the diabatic heating per unit time per unit mass, 

MTS is the forced vertical velocity due to flow over mountain topography, which will 
be discussed in Chapter III, 
f is the Coriolis parameter, 

P is df/dy, 

d() 0() ao ao 5() 

— = — + u — + V— + w — , and 
dt d\. dx dy dZ 

1 is the horizontal Laplacian operator. 

The prognostic equations are linearized by expanding the variables into their 
mean and perturbation states, as in Jordan (1985). The resulting linearized forecast 
equations are: 



— = -f D' - u— - pv' , 
ot dx 



(2.5) 



^D' dD‘ dw du d^(p‘ 

^t dx dx dZ d^ 



( 2 . 6 ) 
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RTw + MTS' , 



(2.7) 



£9s'^ __£9s1_ 

di dx dy 

dT JT _ dT _ dJ 

dt dx 8y dZ ^ 



( 2 . 8 ) 



where R is the gas constant for air, (') denotes perturbation quantities and ( — ) denotes 
mean quantities. The use of Xbar in the text will be used to denote mean quantities of 
a variable X. 

The diagnostic variables, u', v', w' and are calculated from the forecast 
variables using the definitions of divergence, vorticity, the hydrostatic equation and the 
continuity equation. The relationships are given in equations 2.9 through 2.12. 



^u' 

— = D'. 
cx 




dip' 

dx 



RT' . 



5w' 

D' + — = 0 . 
dZ 



(2.9) 

(2.10) 

(2.11) 

(2.12) 



The use of primes to denote perturbation quantities will be discontinued. All quantities 
used in the remainder of the paper will be perturbation quantities unless otherwise 
noted. 

The mean state is assumed to be in hydrostatic and geostrophic balance. The 
term 5Tbar /dy in the first law of thermodynamics can be evaluated by taking d /dy of 
the hydrostatic equation and substituting for 0(pbar /dy from the geostrophic relation, 
5(pbar /dy = — f ubar. Thus, 
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Geostrophic balance of the mean state at the surface implies 



dy 




(2.14) 



The expressions (2.13) and (2.14) are substituted into equations (2.8) and (2.7), 
respectively. 

A singlewave spectral representation is used in the x-direction, with wave number 
p = 271/L, where L is the wavelength in the x-direction. The perturbation quantities 
have the form: 



(^(x,Z,t) = Aj(Z,t) cosux + A 2 (Z,t) sinpx , 



(2.15) 



D(x,Z,t) = Dj(Z,t) cospx + D 2 (Z,t) sinpx , 



(2.16) 



T(x,Z,t) = Tj(Z,t) cospx + T 2 (Z,t) sinpx , 



(2.17) 



(Pg(x,Z,t) = Sj(Z,t) cospx + S 2 (Z,t) sinpx , 



(2.18) 



u(x,Z,t) = Uj(Z.t) cospx + b' 2 (Z,t) sinpx , 



(2.19) 



v(x.Z,t) = Vj(Z,t) cospx + V 2 (Z,t) sinpx , 



( 2 . 20 ) 



w(x,Z,t) = Wj(Z,t) cos|ix + W 2 (Z,t) sin|ix , 



( 2 . 21 ) 



(p(x,Z,t) = H|(Z,t) cos|ix + H 2 (Z,t) siniix , (2.22) 



Q(x,Z,t) = Qi(Z,t) cosp.x + Q 2 (Z,t) sinjix , (2.23) 



MTS(x,Z,t) = MTSj(Z,t) cosjix + MTS 2 (Z,t) sin|ix . 



(2.24) 



The relations (2.15) through (2.24) are substituted into equations (2.5) through 
(2.12). The prognostic and diagnostic equations are separated into equations for the 
cosine and sine terms. The resultant prognostic equations are: 



dt 



-f Dj - upA2 - pVj , 



(2.25) 



dA^ 

dt 



— fD2 + ujiAj — (3V2 , 



(2.26) 



fA, - 1 THD 2 - PU| - , (2.27) 



^2= rAj + -unD, - puj 4 . n ^ W, + , 



(2.28) 



5T, f du 61 

— ^ = -upT^ + V, - — W, + Q, , 



(2.29) 
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(2.30) 



dT') f du ^T 

= upTj + --- + Q2 . 
€t ^ R dZ ^ cZ ^ 


(2.30) 


^S, 

— ^ = -upS2 + f uVj - RTW^ + MTSj . 


(2.31) 


— 

— ^ = upS, + f uV^ - RTW^ + MTS^ . 
dt 1 - ^ ^ 

The resultant diagnostic equations for u and v are: 


(2.32) 


U = - ^2 

' n ■ 


(2.33) 


D| 

1 T _ 1 

, 


(2.34) 


Vj = - — 2 , 


(2.35) 


< 

II 

1^ 


(2.36) 



Geopotential values above the surface are obtained by integrating the hydrostatic 



equation from the surface (Z = Z^) to height Z: 




Z 

Hj = R J T|(Z.t) dZ + Sj . 
Zo 


(2.37) 


Z 

H., = R J T.,(Z,t) dZ + S., . 
Zo “ 


(2.3S) 



The vertical velocity is calculated by integrating the continuity equation from the top 
of the atmosphere (Z = Zj) down to height Z. The upper boundary condition, w = 0 
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at Z = Zj, is used. This boundary condition is not exact, but some form of it is used 
in most numerical models. The diagnostic equations for the vertical velocity are: 

Z 

W| = J Di(Z.t) dZ , 



Z 

W 2 = f D2 (Z,i) dZ . 

Zo 

Equations (2.25) through (2.40) are the prognostic and diagnostic equations that 
govern all four numerical models. Using the given basic state and the one-wave spectral 
perturbation quantities, the governing equations reduce to functions of Z and t. The 
models are effectively one-dimensional. 

To display the results of each model, the sine and cosine amplitudes of each 
variable are combined to determine the amplitude and phase of a single cosine wave in 
the x-direction. A typical variable has the form: 

Y(x,Z,t) = A(Z,t) cos(px-6) , (2.41) 

where the amplitude is A(Z,t) and the phase is 6(Z,t). The amplitude and phase are 
calculated at each level for all variables. 

C. TIME DIFFERENCING 

Two fon^'ard time steps are taken to start each model and then leapfrog time 
differencing is used. The leapfrog scheme is employed because of its ease to code. A 
Robert filter is used to reduce the amplitude of the computational mode generated by 
the leapfrog time differencing. The filter is discussed by Haltiner and Williams (1980). 
For a prognostic variable F, calculate Fbar — the average value of F at time step 
(n— l)At, using equation (2.42), 



(2.39) 

(2.40) 



^n — 1 ^n — 1 *^^^n ^^n — 1 ^n — 2^ ’ 



(2.42) 



where y is a weighting function. Using the unaveraged values at time step nAt, 
compute the tendency from its predictive equation. The predicted value at 

time step (n+ l)At is then calculated using equation (2.43), 
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(2.43) 



^n+l = 

In all the experiments, y = 0.05 is used. The time step for each experiment is 
calculated in the model by requiring, for computational stability, 

1 

vAt = - , (2.44) 

where v = pc, and c is the typical phase speed of an external gravity wave. 

D. VERTICAL GRIDS 

Each of the models uses one of two vertical grids. The two ways used to 
distribute the variables over discrete levels are depicted in Fig. 1.1. The staggered levels 
are represented by the dashed lines in Fig. 1.1. Notice that the heights at which the 
variables are defined change between the two grids. The notation used in this paper to 
denote the staggered and unstaggered levels is consistent with the conventions used in 
the coded models. The height of the unstaggered levels is denoted as Z'. The height of 
the staggered levels is denoted as Z. In the models, both Z'j and Z| are defined to be 
the surface of the earth. It is assumed that the staggered level Z- is exactly in the 
middle of the layer between Z'pj and Z'^. This distinction is important because the 
models can have layers with unequal depth. Thus, the height of the staggered levels is 
defined relative to the height of the unstaggered levels. 

A finite difference model is written for each of the grid structures. The models 
are denoted as FDM-A and FDM-B. Similarly, finite element models using the two 
grids are indicated by FEM-A and FEM-B. 

E. FINITE DIFFERENCE MODELS 

The only differences in the equations between the two FDM models are the 
approximations of terms involving dubar/dZ and ^Tbar/5Z in the prognostic equations 
and the approximations of the integral in the diagnostic geopotential equation. 
Centered difference approximations are used, except at the boundaries where one-sided 
differences are employed. The finite difference approximations used in the prognostic 
equations are listed in Appendix A. 
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F. FINITE ELEMENT MODELS 



1. FEM-A 

The FEM-A model defines vertical velocity (w) at the unstaggered levels in 
terms of the basis functions \|/j(Z). The other variables are defined at the staggered 
levels in terms of the basis functions (Pj(Z). The expansion for a typical term is 

n+ 1 

A,(Z,t) = 1 A^t) vXZ) . (2.45) 

i=i 

The basis functions for this model are depicted in Fig. 2.1. The basis functions V|/(Z) 
are defined for the unstaggered levels (solid lines at height Z') and the basis functions 
(p(Z) are defined for the staggered levels (dashed lines at height Z). 




Fig. 2.1 Basis functions for grids A and B. 
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The finite element approximations for - the vorticity, divergence and 
thermodynamic equations are derived by substituting the expansion for each dependent 
variable into equations (2.25) through (2.30). Each equation is multiplied by (Pj(Z) and 
integrated with respect to Z from the bottom to the top of the atmosphere. Each term 
in the equations is the finite sum of separate integrals. Only the integrals of 
overlapping basis functions are nonzero. The resultant equations, listed in Appendix 
B, are matrix equations. For an N-layer model, the vectors, Aj, A 2 , Dj, D^, H|, H 2 , 
Q], Q 2 . T 2 , U|, L' 2 , V|, and V 2 , contain n+2 components. The vectors Wj. and 

W 2 contain n+ 1 components. The matrices M, K, and <I>, defined below, are 
(n + 2) X (n+ 2) matrices. The matrix P’ defined below, is an (n+ 1) x (n+2) matrix. 

The mass matrix IVl for this model is defined by 



The matrix K is defined for terms multiplied by u, 
i+ 1 Zj 

K,.(u) = y -uk I <p.(Z) (P^(Z) <Pj(Z) dZ for li - j| 1 . (2.47) 




for |i — j| ^ 1 . (2.46) 



k=i-l Z 



0 



du 



dT 



The matrix P is defined for terms multiplied by — W, or by — W, 

dZ dZ 




Vj(Z) (Pj(Z) dZ for|i-i|<l. (2.48) 



where x = u, or T. 



du 



The matrix O is defined for terms multiplied by — V, 



dZ 
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i+ 1 



"T d(p 



d>.j(u) = £ u*^ ^ <p.(Z) <P|(Z) dZ - for |i - j| S 1 . 



k = i-l Z, 



(2.49) 



The staggered basis functions present two general problems for evaluating the 
elements of the four matrices. First, for an n-layer model, portions of basis functions 
(pj(Z) and (Pj^^2(^) defined in the model atmosphere but the physical meaning of 
contributions from those terms is unclear. The contributions are included in the first 
two rows and the last row of each matrix. Second, only portions of basis functions 
<P2(Z) and |(Z) are defined in the model atmosphere. To describe the incomplete 
sides of both basis functions an assumption must be made about the value of (p2 ^he 

surface and j at the top of the atmosphere. 

Assumptions are made and procedures are developed on an attempt to resolve 
these problems. In this model the mean state variables, ubar and Tbar, are defined only 
at the n staggered levels. However, ubar and Tbar values defined at the nodal points 
of (pj(Z) and + important in the Galerkin formulation of the dubar/dZ and 

^Tbar/5Z terms. In these experiments, the values of ubar and Tbar are defined at the 
surface and top of the atmosphere. Jordan (1985) did not define them at the nodal 
points of (pjCZ) and <Pn + 2(^)- of the major modifications of these experiments is 
to define ubar and Tbar at the nodal points of <pj(Z) and 4-2(2^)- constant 
shear with height, ubar and Tbar are defmed at the boundaries such that the shear in 
the two half layers at the boundaries is the same as the shear in the other layers. To 
evaluate the staggered basis functions defined in the layers between the surface and Z2, 
and Zj^^j and the top of the atmosphere, it is assumed that the value at the 
boundaries of those basis functions is one-half Thus, three-fourths of the basis 
functions (P2(^) ^n + 2^^^ defmed in the model atmosphere. 

The equations for the general elements of the four matrices are evaluated by 
substituting into equations (2.46) through (2.49) the formulas for (p- ^.^(Z), (p-(Z), 
(Pj_2(Z), \J/^^j(Z), Vj(Z), \}/-_j(Z), and \j/^_2(Z), in terms of the local coordinate 
^ = Z — Z-. The equations for these basis functions defined for the levels 1, 2, i, and 
n+ 1, are listed in Appendix C. The matrices were evaluated by integrating numerically 
using 2 point Gausian Quadratures. 

The vorticity, divergence and thermodynamic equations, written in matrix and 
vector form, are: 
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(2.50) 



dA, 

M^l= M(-fDi - p Vj) - ^iK(u)A2. ' 

dA-) 

M — M(-fD 2 - p V^) + H K(u) Aj , 
dD I -) 

M — ‘= M(fA, - p U, + K(u) H P(u) W, , 

dt 111 L 1 

dD-) -) 

M — M(f A 2 - p U 2 + |i^H 2 ) + [i K(u) Dj - ji P(u) W 2 , 
dT, f 

M — * = - K(u) + - 1>(u)Vj - P(T) W, + MQ, , 

dt R 

dT r “ 

M — K(u) Tj + -<I>(u)V 2 - P(T) W 2 + MQ 2 . 



(2.51) 

(2.52) 

(2.53) 

(2.54) 

(2.55) 



Equations (2.50) through (2.55) are simplified by multiplying each equation by 
M*^ and applying the Robert filter. Actually one should not compute the inverse of 
M. Instead at t = 0 one should obtain the LU factorization of M. Thus at each time 
step one only needs to forward and back solve a triangular system. The matrices 
M*^ K, M-^ P, and M'^ <t> are constants. They are constructed in the initialization 
subroutine and stored for use in the forecast subroutine. The matrices are multiplied 
by the appropriate vectors with values for time level nAt. The resultant forecast 
equations are vector equations and the forecast value for the i-th vertical level is the 
sum of values in the i-th location of each vector equation. The prognostic equations 
for the vorticity, divergence and potential temperature vectors are: 



Al(n+I) = + 24t(-fD[ - pV, - hM-'K(u)A 2 )(„) , (2.56) 

A2(n+1) = + 2M(-fD2 - pV^ + hM-'K(u)A,)(„) , (2.57) 
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(2.58) 



^l(n+ 1) “ ^l(n- 1) 

+ 2At(fAj-pUj + - ^iM*^K(u)D 2 - ^lM■^PCu)W2)(n) 

^2(n+l) " ^2(n-l) 

+ 2At(fA2-pU2 + u^H 2 - ^lM*^K(u)D| + iiiVr^P(u)W|)^j^^ , (2.59) 

"^Kn+l) "" ^l(n-l) 

+ 2At (M'^K(u)T 2 + -M-^<D(u)Vj - M'^P(T)Wj + Q\\r{) - (2-60) 

^2(n+l) " '^2(n- 1) 

+ 2At (-M'^K(u)Tj +^M'^(D(u)V 2 - M'^P(T)W 2 + Q2)(n) > (2.61) 

where the subscripts (n+ 1), (n) and (n — 1) refer to the values of the vectors at time 
step (n+ l)At, nAt and (n— l)At, respectively. 

The surface geopotential and the diagnostic variables are calculated using the 
corresponding equations in model FDM-A (see Appendix A). 

2. FEM-B 

The FEM-B model defines vertical velocity, potential temperature, mean state 
potential temperature and diabatic heating at the unstaggered levels in terms of the 
basis functions \|/j(Z). The other variables are defined at the staggered levels in terms 
of the basis functions <Pj(Z). The basis functions are the same as defined for the FEM- 
A model, shown in Fig. 2.1. 

The finite element approximations for the vorticity, divergence and 
thermodynamic equations are derived by substituting the expansion for each dependent 
variable into equations (2.25) through (2.30). The vorticity and divergence equations 
are multiplied by (Pj(Z) and integrated with respect to Z from the bottom to the top of 
the atmosphere. The resultant Galerkin formulation of the vorticity and divergence 
equations are the same as those derived for model FEM-A. The matrices in those 
equations, M, K, and P, are given by equations (2.46) through (2.48) as defined for 
FEM-A. The thermodynamic equations are multiplied by M/j(Z) because potential 
temperature is defined at the unstaggered levels. As before, the equations are integrated 
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through the depth of the atmosphere. The resultaat- equations are listed in Appendix 
D. Four additional matrices are defined for the two thermodynamic equations. The 
mass matrix FI is 




Z 



V|/.(Z) Vi(Z) dZ 



for|i-j|^l. (2.62) 



0 



The matrix F is defined for terms multiplied by u, 



r„(u)=y;u''j V.(Z) <p^(Z) Vi(Z) dZ for|i-j|<l, (2.63) 

k = i-l Z„ 



i+ 1 Zy 



Terms multiplied by — V, give rise to the transpose of the matrix P, defined by (2.48). 



As discussed in the FEM-A model description, the staggered finite elements 
present problems for evaluating the elements of the matrices. In this model ubar is 
defined at the surface, the top of the atmosphere, and at the n staggered levels. The 
mean state temperature, Tbar, is defined at the unstaggered levels so special definitions 
for it are not needed. Jordan (1985) did not include the contributions from the 
perturbation quantities defined at the nodal points of (Pj(Z) and + 
included in this model as pan of the modifications of this experiment. The staggered 
basis functions, (pj(Z), are evaluated at the boundaries using the assumptions discussed 
in the previous section. 

The elements of matrices 11, F, and T are evaluated by substituting formulas 



dT 



The matrix T is defined for terms multiplied by — W, 



dZ 
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terms of the local coordinate ^ = Z — Z'j, into equations (2.62), (2.63), and (2.65). 
Formulas for these basis functions are listed in Appendix C. As in FEM-A the 
matrices are evaluated by integrating numerically using 2 point Gausian Quadratures. 

The forecast matrix equations for vorticity, divergence and temperature are 
simplified in a manner similar to the method described for model FEM-A. The final 
form of the vorticity and divergence vector equations are the same as for model FEM- 
A, equations (2.56) through (2.59). The thermodynamic vector equations are: 

^l(n+ 1) ” ^l(n- 1) 

+ 24t ( + nn-‘r(u)T2 + -^n->pT(u)V, - n->4'(f)W, + , (2.66) 

^2(n+ 1) “ ^2(n- 1) 

+ 26t (-nn-‘r(u)T, + f n-‘pT(u)V 2 - n-‘'P(f)W 2 + Q 2 )(„) . (2.67) 

In this model the vectors, Aj, A 2 , Dj, D 2 , Hj, H 2 , Up U 2 , Vp and V 2 , 
contain n+2 components. The vectors Qp Q 2 , Tp T 2 , Wp and W 2 contain n+ 1 
components. The matrices IT, F, and T are (n+ 1) x (n+ 1) matrices and the matrix 
is an (n+2) x (n+1) matrix. The surface geopotential and the diagnostic 
variables are calculated using the corresponding equations in model FDM-B. 
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III. EXPERIMENTS AND RESULTS 



Four experiments are performed with each model; an initial perturbation in the 
meridional flow, flow over mountain topography, flow with a diabatic heat source and 
baroclinic flow with vertical shear in the ubar field. The first two and the last 
experiments are performed with six- and then sixty-layer models. The heating 
experiment is repeated with six-, twelve- and sixty-layer models. The analytic solution 
of each experiment has not been derived. For each experiment, the sixty-layer model 
results are intercompared to determine if the models are converging to the same 
solution. The standard of comparison for each six-layer model is its corresponding 
sixty-layer solution. Temperature and divergence profiles are examined in each 
experiment. 

Several parameters are defined identically in each experiment. The vertical 
coordinate, Z is defined between zero and one (1000— 368mb) and the vertical levels are 
equally spaced. The x-wavelength is 4,000 kilometers. The time step is 17.7 minutes. 
The Coriolis parameter is defined at 45 degrees latitude. The mean state potential 
temperature increases with height from its surface value of 310.0 degrees K (Kelvin). 

A. ROSSBY WAVE EXPERIMENT 

Rossby waves are generated in each model using an initial perturbation, v' = 5.0 
meters/ second (m/s), in the cosine term of the meridional flow. All other perturbations 
are initially zero. There is no diabatic heat source and no mountain topography. There 
is no vertical shear in the ubar field and ubar = 10.0 m/s. The latitudinal variation of 
the Coriolis parameter, P, is defined at 45 degrees latitude. The forecast experiments 
are terminated at 96 hours. 

1. Sixty- Layer Models 

The sixty-layer FDM-A, FDM-B, and FEM-B models converge to the same 
temperature and divergence solutions (Figs. 3.1— 3.2). All figures will be found at the 
end of the chapter. Note that the staggered models FDM-A and FEM-A are defined 
as zero at their lowest level because this level is' below the bottom of the atmosphere. 
It should also be noted that the phase of each variable is defined between zero and 360 
degrees. There is a discontinuity in the phase profile if the phase passes through zero 
degrees. The height at which the temperature phase discontinuity in model FDM-A 



28 



occurs differs from the other two models because temperature is defined at the 
staggered levels in FDM-A. The three models represent the same physical solution, 
which is called the consensus solution. 

The FEM-A temperature amplitude is slightly smaller than the consensus 
amplitude, 2% at height Z = 0.10 , and an amplitude oscillation is present in the 
lowest three layers of the atmosphere. In the results of Jordan (1985) (Fig. 3.3(top)), 
the FEM-A solution had a 32% difference from the consensus solution at height 
Z = 0.10 and a much more jagged profile in the lowest three layers of the atmosphere. 
Jordan (1985) felt that this jagged profile may be caused by the terms in the matrices 
which represent contributions from the basis functions near the low^er boundary of the 
model. The FEM-A model is modified near the boundary to correct this problem. The 
FEM-A temperature phase is within 0.1 degree of the consensus. It is high enough to 
pass through zero one level above the consensus. The shape of the FEM-A divergence 
amplitude is consistent with that of the consensus amplitude, however it is slightly 
lower at the bottom and slightly higher at the top of the atmosphere. The consensus 
divergence phase is nearly constant with height and the divergence phase profile for 
model FEM-A is almost identical with the consensus profile. 

In the results of Jordan (1985) (Fig. 3.3(bottom)) the FEM-B solution had a 
5% difference from the consensus solution and a major oscillation in the lower layers 
of the atmosphere. The FEM-B model is modified near the boundary to correct this 
problem and is now part of the consensus solution (Figs. 3.1— 3.2) 

2. Six-Layer Models 

The comparison of six- and sixty-layer profiles for variables defined at 
staggered levels may be initially misleading. The first staggered level in a six-layer 
model occurs at Z = 0.0833. The lowest staggered level in a sixty-layer model is 
defined at Z = 0.0083, which may be mistaken for the surface in the graphs. When the 
values of a six-layer model coincide with a sixty-layer profile, the models are considered 
to represent the same physical solution, even though the six-layer model has a smaller 
vertical domain for staggered variables. 

The temperature amplitude profiles of the six-layer models are very similar to 
their corresponding sixty- layer results (Figs. 3.4— 3.5). The previously discussed 
problems in the lowest three layers of model FEM-A are still slightly evident in the 
lowest three layers of the six-layer profile (Fig. 3.5(top)). The six-layer profile of the 
FEM-A model is nevertheless a very good approximation of the consensus profile. 
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The six-layer divergence amplitude profiles for the grid A models approximate 
their sixty-layer counterparts in a similar manner (Figs. 3.6(top)— 3.7(top)). And the 
six-layer divergence amplitude profiles for the grid B models approximate their sixty- 
layer counterparts in a similar manner also (Figs. 3.6(bottom)— 3.7(bottom)). The grid 
B models more closely approximate the sixty-layer consensus profile in the lower half 
of the atmosphere and the grid A models more closely approximate the sixty-layer 
consensus profile in the upper half of the atmosphere. The six-layer profile of the 
FEM-A model is the closest of the four in their approximation of their sixty-layer 
counterparts and the sixty-layer consensus profile (Fig. 3.7(top)). 

B. MOUNTAIN TOPOGRAPHY EXPERIMENT 

The forced vertical velocity term, MTS, in the surface geopotential forecast 
equation (2.3) is non-zero in this experiment. It represents the contribution to surface 
geopotential from air flowing over mountain topography which varies sinusoidally in 
the x-direction and has no variation in the y-direction. The mountain ridge-to-valley 
height difference is 1,500 meters. To reduce the trauma for the model, the mountains 
are gradually "built" to their full height over a period of 36 hours. Thus, the forced 
vertical velocity increases in the first 36 hours of the forecast period and is constant for 
the remainder of the 96-hour forecast period. The equations used to define the forced 
vertical velocity are included in Appendix E. There is no vertical shear in the ubar field 
and ubar = 10.0 m/s. (3 and all other initial perturbations are zero in this experiment. 

1. Sixty- Layer Models 

The FDM-A, FDM-B, and FEM-B models converge to the same physical 
solution for temperature and divergence (Figs. 3.8 and 3.9). 

The FEM-A temperature amplitude is again slightly smaller than the 
consensus amplitude (2%) and a small oscillation is again present in the lowest three 
layers of the atmosphere. The FEM-A temperature phase and divergence phase 
solutions are exactly the same as the consensus solutions. The FEM-A divergence 
amplitude is only very slightly lower than the consensus amplitude at the bottom of the 
atmosphere (0.5%) and is indistinguishable from the consensus profile by Z = 0.20. 
Jordan (1985) found that the temperature amplitude of the unmodified FEM-A model 
was 30% higher than the consensus near the bottom of the atmosphere (Fig. 
3.10(top)) and that the divergence amplitude of the unmodified FEM-A model was 
higher at the bottom of the atmosphere and lower at the top of the atmosphere than 
the consensus (Fig. 3.10(bottom)). 
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Jordan (1985) again found a jagged temperature amplitude profile in the 
unmodified FEM-B model in the lowest two layers of the atmosphere (Fig. 3.11(top)) 
and the temperature amplitude w’as 5% more than the consensus near the surface. The 
unmodified FEM-B temperature phase profile was jagged in the lowest three layers and 
top two layers of the atmosphere while the rest of the profile was within one degree of 
the consensus (Fig. 3.1 l(bottom)). The modified FEM-B profiles are now part of the 
consensus in every case for this experiment (Figs. 3.8 and 3.9). 

2. Six-Layer Models 

The temperature amplitude profiles of the six-layer models FDM-A, FDM-B, 
and FEM-B are identical with each other and also with the consensus solution (Figs. 
3.12(top) and 3.13(bottom)). The slight oscillation in the lowest three layers of the 
FEM-A model is still slightly evident in the six-layer model (Fig. 3.13(top)), however 
the six-layer model profile is still very close to the consensus profile (Fig. 3.12(top)). 

The six-layer divergence amplitude profiles are very close to the sixty-layer 
consensus (Figs. 3.14 and 3.12(bottom)). Both of the six-layer finite element model 
divergence amplitude profiles are closer to the sixty-layer consensus divergence 
amplitude profile than the FDM-A profile (Figs. 3.14 and 3.15) 

C. DIABATIC HEATING EXPERIMENT 

A diabatic heat source is defined in the layer between Z = 0.40 and Z = 0.60 
(670— 549mb). The rate of heating is constant in time and varies in x and Z, 

z — z 

Q(x,Z,t) = HEATING cos^( ^ n ) cos()ix) , (3.1) 

where HEATING is 5.0 K/day, is the midpoint of the heated layer, and Z^^ and Z^ 
are the low'er and upper boundaries of the heated layer, respectively. The diabatic 
heating vectors, Qj and Q 2 , are defined in the initialization subroutine and stored for 
use in the forecast subroutine. There is no vertical shear in the ubar field and ubar = 
10.0 m/s. P and all other initial perturbations are zero in this experiment. The forecast 
length is 96 hours, however, a 12-hour forecast is made for comparison with Jordan 
(1985). 

1. Sixty- Layer Models 

For the diabatic heating function defined in equation (3.1), the maximum 
heating occurs at Z = 0.50, the midpoint of the heated layer. The models defined 
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using grid B define temperature and the heating functions at this point. The grid A 
models do not have temperature and diabatic heating defined at this point so the 
maximum rate of heating in these models is slightly less than in the other models, and 
the maximum heating occurs throughout one layer rather than occuring at one point. 
The heating rate at each level is listed in Appendix F for the six-, twelve- and sixty- 
laver models. 

The sixty-layer profiles for the four models are quite similar, the differences 
occur mainly because the models are responding to different forcing. The temperature 
amplitude profiles for the B grids come to a sharp point at Z = 0.50 and the grid A 
models have a square-nosed profile around this point (Fig. 3.16(top)). The previously 
identified temperature amplitude oscillations in the lowest layers of models FEM-.A and 
FEM-B were not evident in the unmodified models (Jordan, 1985) (Fig. 3.17) and are 
not evident in any of the modified models (Figs. 3.16(top)). This is because the heating 
is defined far enough away from the boundaries so that the previous problems with the 
basis functions at the bottom and the top of the atmosphere do not show up here. 
Note that there is a jagged profile at the bottom and top of the heated layer in the 
twelve hour forecasts because they are not steady state yet. In summary, the sixty- 
layer temperature profiles of all four models represent the same physical response to 
the diabatic heating. 

The shape of the divergence amplitude profile is not symmetric around 
Z = 0.50 because divergence is not defined exactly at the midpoint of the heated layer 
in either of grids A or B (Fig. 3.18(top)). The divergence amplitude is identical outside 
the heated layer for all models except FEM-A which was also true at 12 hours with the 
unmodified models (Jordan, 1985), (Fig. 3.19(top)). The divergence phase profiles are 
virtually identical for all models in this experiment. Jordan (1985) found that the 
divergence phase for the unmodified FEM-A model at 12 hours was slightly different 
from the consensus outside the heated layer (Fig. 3.19(bottom)). The modified FEM-A 
model profile at 12 and 96 hours is exactly the same as the consensus profile (Figs. 3.20 
and 3.1S(bottom)). 

2. Six and Twelve-Layer Models 

The difference between grids is more evident in this experiment than in the 
other experiments. Each model is run with both six and twelve layers and the results 
are compared with the sixty-layer consensus of the four models. 
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The six-layer grid A model temperature and divergence amplitude fields barely 
respond to the diabatic heating (Figs, 3.21(top) and 3.22(top)). The grid A models are 
only about 6% of the amplitude of the grid B models. Comparing the maximum 
diabatic heating terms for the two grids, in Appendix F, it is found that the maximum 
heating in the grid A models is only about 6% of the maximum heating in the grid B 
models. This is because heating is defined to be greatest at the midpoint of the heated 
layer and decrease as you get away from the midpoint. Grid B defines heating at the 
maximum, or midpoint and grid A defines heating at points equidistant from the 
midpoint. As the number of layers gets smaller the distance between the midpoint and 
the first layer in grid A which is used to define heating increases, thus decreasing the 
value of the heating that goes into the layer. 

The grid A models have a stronger response to heating with twelve layers than 
with six layers because with twelve layers heating is defined closer to the maximum at 
the midpoint of the heated layer (Figs. 3.21—3.22), The six- and twelve-layer models 
converge to the sixty-layer profile equally well for the finite element models as for the 
finite difference model for grid B (Figs. 3.23 and 3,24). For grid A the twelve-layer 
divergence amplitude is closer to the corresponding 60 layer profile for the finite 
element than for the finite difference model (Fig. 3.25) but the finite difference twelve- 
layer profile is closer to the consensus (Fig. 3.26). 

D. BAROCLINIC INSTABILITY EXPERIMENT 

Vertical shear in the ubar field is defined in this experiment. The wind profile is a 
linear function of Z, with ubar = (STRGTH) Z, where STRGTH defines the strength 
of the wind at the top of the atmosphere in m/s. In this experiment STRGTH is 
defined as 40 m/s. Waves are generated in each model using an initial perturbation, 
v' = 5.0 m/s in the cosine term of the meridional flow. P and all other perturbations 
are initially zero. There is no diabatic heat source and no mountain topography. The 
forecast experiments are terminated at 96 hours. 

1. Sixty- Layer Models 

The sixty-layer FDM-A, and FDM-B models converge to the same solution 
for the temperature amplitude, temperature phase, divergence amplitude and divergence 
phase (Figs. 3.27 and 3.28), 

The FEM-A model shows a slightly higher temperature amplitude than the 
consensus profile but the shape is the same as the consensus profile. The same 
problem is evident for the divergence amplitude profile but the magnitude of the error 
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is smaller. The FEM-A model shows no sign of the oscillation that is evident near the 
surface in the Rossby wave experiment and the mountain topography experiment. 
This may be because the ubar field is almost zero near the bottom of the atmosphere 
where Z is small which cancels out the problem with the basis function that may be 
causing the oscillation in the other experiments. The FEM-A model profile is almost 
identical with the consensus for temperature phase and identical with the consensus 
profile for divergence phase. 

The FEM-B model profile shows a lower temperature amplitude than the 
consensus profile and the shape is the same except for a small oscillation that occurs in 
the lowest three layers of the atmosphere. This oscillation is not evident in any of the 
other experiments. For divergence amplitude, the amplitude is lower than the 
consensus profile, but there is no oscillation in the lower layers of the atmosphere and 
the magnitude of the error is smaller, as it is for the FEM-A model. The FEM-B 
model profile is further from the consensus profile than the FEM-A model profile for 
temperature phase but it is still very close (Fig. 3.29). The profiles for divergence 
phase are identical for all four models (Fig. 3.28(bottom)). 

2. SLx-Layer Models 

In the six-layer case the FDM-A, FDM-B, and FEM-A models all seem to 
converge to the same solution, but the FEM-B profile is off in the temperature 
amplitude and the divergence amplitude (Figs. 3.30 and 3.31) although the pattern is 
similar. The phase profiles for all four six-layer models are very close to the sixty-layer 
consensus. The FEM-A six-layer model profile for temperature amplitude is closer to 
the sixty-layer consensus than the FDM-A six-layer model and slightly better than the 
FD.M-B six-layer model (Fig. 3.32). The temperature amplitude profile oscillation in 
the lowest layers of the FEM-B sixty-layer model is evident at Z = 0.3 in the six-layer 
model (Fig. 3.33). The reason for this oscillation is unclear, but it should have 
something to do with the basis functions near the boundaries. 
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TEMPERATURE AMPUTUDE (KELVIN) 



TEMP. PHASE V = 5.0 CASE (60-L) 



o 




Fig. 3.1 Sixty-layer Rossby wave experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.2 Sixty-layer Rossby wave experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) arc compared. 
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Fig. 3.3 Sixty-layer Rossby wave experiment at 96 hours from Jordan (1985). 
Temperature amplitude profiles are compared for models FEM-A (top) and 
FEM-B (bottom) and FD.M-C, which represents the consensus profile. 
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TEMP. AMP. V=5 CASE (PDM-A) 
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Fig. 3.4 Six-layer Rossby wave experiment at 96 hours. 
Temperature amplitude profiles are compared for the six-layer 
and sixty-layer FD.M-A (top) and FD.M-B (bottom) models. 
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Fig. 3.5 Six-layer Rossby wave experiment at 96 hours. 
Temperature amplitude profiles are compared for the six-layer 
and sixty-layer FEM-A (top) and FEM-B (bottom) models. 
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Fig. 3.6 Six-layer Rossby wave experiment at 96 hours. 
Divergence amplitude profiles are compared for the six-layer 
and sixty-layer FDM-A (top) and FDM-B (bottom) models. 
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Fig. 3.7 Six-layer Rossby wave experiment at 96 hours. 
Divergence amplitude profiles are compared for the six-layer 
and sixty-layer FE.M-A (top) and FE.M-B (bottom) models. 
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Fig. 3.8 Sixty-layer mountain topography experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.9 Sixty-layer mountain topography experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) are compared. 
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Fig. 3.10 Sixty-layer mountain topography experiment at 96 hours from Jordan (1985). 
Temp, (top) and divergence (bottom) amplitude profiles are compared 
for models FEM-A and FDM-C, which represents the consensus profile. 
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TEMPERAIURC phase FOR MOUNTAIN CASE (60 UYERS) 




Fig. 3.1 1 Sixty-layer mountain topography experiment at 96 hours from Jordan (1985). 
Temperature amplitude (top) and phase (bottom) profiles are compared 
for models FEM-B and FDM-C, which represents the consensus profile. 
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Fig. 3.12 Six-layer mountain topography experiment at 96 hours. 
Temperature amplitude profiles (top) and divergence amplitude 
profiles (bottom) are compared. 
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TEMP. AMP. MOUNTAIN CASE (FEM-A) 
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Fig. 3,13 Six-layer mountain topography experiment at 96 hours. 
Temperature amplitude profiles are compared for the six-layer 
and sixty-layer FEM-A (top) and FEM-B (bottom) models. 
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Fig. 3.14 Six-layer mountain topography experiment at 96 hours. 
Divergence amplitude profiles are compared for the six-layer 
and sixty-layer FE.M-A (top) and FE.M-B (bottom) models. 
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Fig. 3.15 Six-layer mouniain topography experiment at 96 hours. 
Divergence amplitude profiles are compared for the six-layer 
and sixty-layer FD.M-A model. 
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TEMP. AMP. HEATING CASE (60-L) 
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Fig. 3.16 Sixty-layer diabatic heating experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.17 SL\ty-layer diabatic heating experiment at 12 hours from Jordan (1985). 
Temperature amplitude profiles are compared for models FEM-A (top) 
and FEM-B (bottom) and FDM-C, which represents the consensus profile. 
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Fig. 3.18 SLxty-laycr diabaiic heating experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) are compared. 
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Fig. 3.19 Sixty-layer diabatic heating experiment at 12 hours from Jordan (1985). 
Divergence amplitude (top) and phase (bottom) profiles are compared 
for models FEM-A and FDM-C, which represents the consensus profile. 
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Fig. 3.20 Sixty-layer diabatic heating experiment at 12 hours, 
divergence phase profiles are compared for models FEM-A 
and FDM-A, which represents the consensus profile. 
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TEMP. PHASE HEATING CASE (6-L) 
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Fig. 3,21 Six-layer diabatic heating experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.22 Six-layer diabaiic healing experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) are compared. 
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TEMP. PHASE HEATING CASE (12-U 
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Fig. 3.23 Twelve-layer diabatic heating experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.24 Twelve-layer diabatic heating experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) are compared. 
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Fig. 3.25 Twelve-layer diabatic heating experiment at 96 hours. 
Divergence amplitude profiles are compared for the twelve-layer 
and sixty-layer FDM-A (top) and FE.M-A (bottom) models. 
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Fig. 3.26 Twelve-layer diabaiic heating experiment at 96 hours. 
Divergence amplitude profiles are compared for models FDM-A, FDM-B, FEM-A, 
FEM-B and sixty-layer FEM-B, which represents the consensus profile. 
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TEMP. AMP. BAROCLINIC CASE (60-L) 
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Fig. 3.27 Sixty-layer baroclinic instability experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.28 Sixty-layer baroclinic instability experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) are compared. 
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TEMP. PHASE BAROCLINIC CASE (60-L) 
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Fig. 3.29 Sixty-layer baroclinic instability experiment at 96 hours. 
Temperature phase profiles are compared on a closer scale. 
Compare Fig. 3.27(bottom). 
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TEMP. PHASE BAROCLINIC CASE (6-L) 
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Fig. 3.30 Six-layer baroclinic instability experiment at 96 hours. 
Temperature amplitude profiles (top) and temperature phase 
profiles (bottom) are compared. 
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Fig. 3.31 Six-layer baroclinic instability experiment at 96 hours. 
Divergence amplitude profiles (top) and divergence phase 
profiles (bottom) are compared. 
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Fig. 3.32 Six-layer baroclinic instability experiment at 96 hours. 
Temperature amplitude profiles are compared for models FDM-A, FD.M-B, FEM-A, 
and sixty-layer FDM-A, which represents the consensus profile. 
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TEMPERATURE AMPUTUDE (KELVIN) 



Fig. 3.33 Six-layer baroclinic instability experiment at 96 hours. 
Temperature amplitude profiles are compared for the six-layer 
and sixty- layer FEM-B model. 
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IV. CONCLUSIONS 



The modification of the finite element models with the inclusion of the basis 
functions near the boundaries ((pj(Z) and (p^+ 2 (Z)) completely fixed the FEM-B 
model in the experiments that can be compared with the results of Jordan (1985). The 
modified FEM-A model is much better when compared with the results of the 
unmodified model, how'ever the unusual temperature amplitude behavior in the lowest 
layers of the model are not totally gone after the modifications. The theory of Jordan 
(1985) that the oscillation in the temperature profiles of both unmodified models is 
generated by matrix elements which represent the contributions from the basis 
functions near the surface, is well supported by the results of the modified models. 
There is still some question as to how to define terms at or below the bottom of the 
atmosphere and a change in the definition of those terms may fix the small oscillations 
that are evident in the FE.M-A model in the Rossby wave and mountain topography 
experiment. 

The finite element models display a better convergence to the sixty-layer 
consensus than the finite difference models in many of the cases and in the other cases 
they are the same. 

Jagged temperature profiles are not observed in the diabatic heating experiment 
in either the modified or the unmodified (Jordan, 1985) models. This may be because 
the heating is defined far enough away from the boundaries that the previous problems 
associated with the boundary terms are not significant. The FEM-A model has a 
slightly different sixty-layer divergence amplitude response outside the heated layer 
than the other models. The differences between grids is most apparent in this 
experiment. The differences may be caused by the difference in the maximum amplitude 
of the heating defined for the different grids. The results may be much closer if the 
maximum heating were to be set equal for both grids. This is supported by the fact 
that the sixty-layer profiles of all four models represent the same physical solution. 

Both of the finite element model results differ from the consensus results for the 
baroclinic instability experiment for temperature and divergence amplitudes. The FEM- 
A model difference is probably a manifestation of the problems observed in the Rossby 
wave experiment and the mountain topography experiment. There is an oscillation and 
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a difference from the consensus in the baroclinic instability experiment for the FEM-B 
model which may also be caused by the definition of terms at or below the bottom of 
the atmosphere. In future studies the temperature on the boundary in the FEM-B 
model should be included in the forecast field. This may improve the behavior for the 
baroclinic experiments, where the surface temperature is very important. 



69 



APPENDIX A 

FINITE DIFFERENCE APPROXIMATIONS 



du 

1. For terms of the form — W : 

dZ 

a. FDM-A and FDM-B, at level Z = Z. : 

’ I 



du 



— W = - (W. (- 

o'- 1 ^ 
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— U: 



dZ 



^i+1 
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-) + Wi_,( ■ _ ) ) 
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2. For terms of the form — W : 

dZ 

a. FDM-A at level Z = Zj : 
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— W = - (W: (- 

A~7 ^ 1 '• 
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b. FDM-B at level Z = 
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du 

3. For terms of the form — V ; 

dZ 

a. FDM-A, at level Z = Zj ; 
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^i ~ ^i-1 

Zi - Zj _ j 



•)} 



b. FDM-B, at level Z = Z'j : 
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APPENDIX B 

GALERKIN FORM OF FEM-A PROGNOSTIC EQUATIONS 



1. Vorticity Equations (2.25) and (2.26) : 




(B.l) 



i+ 1 



dAJ^ 



j = i - 1 dt 



I -f I 



i+ 1 






j = i-l 




i+ 1 

+ [I X 
k=i- 1 





i+ 1 

I vi. 

j=i-l 



<Pj <Pi dZ . 



(B.2) 



Note that in these equations, and the equations that follow, the basis functions 
are functions of Z (<Pj = 9j(Z) and Vj/- = Vj/-(Z)). All of the other variables. A, D, H, 
Q, T, u, U, V, and W, are functions of time ( A- = A-(t), D- = Dj(t), H- = H-(t), 
Qi = Qi(t). T- = Tj(t), u = u(t), Uj = Uj(t), Vj = V.(t), W- = W^(t) ). 
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2. Divergence Equations (2.27) and (2.28) : 






i+ 1 



S 17 Sy = f S A), I '<(>j(Pi 

j = i-l“' ^0 j = i-l ^0 
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(B.4) 
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3. Thermodynamic Equations (2.29) and (2.30) : 

dli, ,Zj 

Z -'iJviVidz = 

j = i_ldt Zq 



i+1 i+1 



- n z z T>2 <Pk <(>i 
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APPENDIX C 

BASIS FUNCTION EQUATIONS FOR FEM-A 



1. Notation: 



Aj Zj Zj _ I 

= z-i -z-i_, 

£ = Z - Z; 



2. For the general case; 



a - A': + .5A':_ 1 
(p.(^) = 1 1 — L 

.5(A'i + 



-A'2 - ^ ^ ^ ^ 

^22 
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2 2 
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3. For special cases, <Pj, (P2, (p^+ 9n + 2» '^n+ I 

-t + .5AN 
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APPENDIX D 

GALERKIN FORM OF FEM-B PROGNOSTIC EQUATIONS 



1. The vorticity equations have the same form as the vorticity equations for model 
FEM-A, equations (B.l) and (B.2). 

2. The divergence equations are the same as the divergence equations in model FEM- 
A, equations (B.3) and (B.4). 

3. Thermodynamic Equations (2.29) and (2.30) ; 




-n I f jij F'’'vj<PkVi'iz 



i+l i+1 



k = i- 1 j = i— 1 “^0 





i = i - 1 ^0 



i+ 1 



(D.l) 
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dTU 
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k= i - 1 j = i- 1 ^0 
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i f^T ^^k 



+ _ y u*' y vJ. f ^ 
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(pj M/| dZ 



-‘fVy Wij J|T J|,k y y. dZ 
k = i-l j = i-l ^0 



i+ 1 



+, I .Q'a Jz’^VjVidz. 



J = l- 1 



(D.2) 



Note that in these equations the basis functions are functions of Z (tp- = <Pj(Z) 
and \|/j = V}/-(Z)). All of the other variables, Q, T, u, V, and W, are functions of time 
(Qi = Qi(t). y = Ti(t), u = u(t), V. = Vi(t), Wi = Wi(t) ). 
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APPENDIX E 

FORCED VERTICAL VELOCITY 



1. The contribution to the surface geopotential from the forced vertical velocity is (Pj, 

(pg(x,t) = (p^ sin^(7tt/2T) sin (px) t < T 

= <Pj^sin(px) t>T, (E.l) 

where (p^^^ is mountain geopotential (m^/s^), t is time, and T is the total time to build 
the mountain. <p^ is a constant, 

^m “ §^m ’ (E-2) 

where g is gravity and is the height of the mountain. is a parameter specified 
in each model. = 750 meters in the thesis experiments. 

2. The time rate of change of <Pj is separated into sine and cosine components for use 
in the surface geopotential forecast equations. 



= .MTSj(t) cos (px) + MTS 2 (t) sin (px) , (E.3) 

d d _ d 
where — = -r- + u.r_ — . 
dt dt dx 



Equation (E.l) is substituted into equation (E.3) and the resultant expression is 
separated into sine and cosine equations. The equations to calculate the terms MTSj 
and MTS 2 are 



MTS|(t) = <(>„, sin^(nt/2T) t^T 

“ Sfc ‘Pm 



and 
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sin (7U/2T) cos (7it/2T) t ^ T 



MTS2(t) 



T 

0 



t > T . (E.5) 



These terms are calculated for each time step in the model's forecast subroutine. 
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APPENDIX F 

DIABATIC HEATING TERMS 



For the diabatic heating function defined in equation (3.1), the maximum heating 
occurs at Z = 0.50, the midpoint of the heated layer. Temperature and diabatic 
heating are defined at the staggered levels for grid A, and at the unstaggered levels for 
grid B. Consequently, the rate of heating differs between the staggered and 
unstaggered models. In these experiments, the heated layer is between Z = 0.40 and 
Z = 0.60, the heating rate is 5.0 K/day, and only the cosine term, Q^, is nonzero in the 
heated layer. The value of the heating term is listed below for staggered and 
unstaggered levels for six-, twelve- and sixty-layer models. 



Grid A 


Six- Layer Models 


Grid B 


Six- Layer Models 


Z 


Qi 


Z 


Qi 


0.250 


O.OOOOOOE + OO 


0.333 


O.OOOOOOE + OO 


0.417 


0.387657E-05 


0.500 


0.578704E-04 


0.583 


0.387668E-05 


0.667 


O.OOOOOOE + OO 


0.750 


O.OOOOOOE + OO 






Grid A 


Twelve- Layer Models 


Grid B 


Twelve- Layer Models 


Z 


Qi 


Z 


Qi 


0.375 


O.OOOOOOE + OO 


0.333 


O.OOOOOOE + OO 


0.458 


0. 36424 1E-04 


0.417 


0.387657E-05 


0.542 


0.364243E-04 


0.500 


0.578704E-04 


0.625 


O.OOOOOOE + OO 


0.583 


0.387668E-05 






0.667 


O.OOOOOOE + OO 


Grid A 


Sixty-Layer Models 


Grid B 


Sixty- Layer Models 


Z 


Qi 


Z 


Qi 


0.392 


O.OOOOOOE + OO 


0.400 


O.OOOOOOE + OO 


0.408 


0.985885E-06 


0.417 


0.387660E-05 


0.425 


0.847474E-05 


0.433 


0.144676E-04 


0.442 


0.214459E-04 


0.450 


0.289351E-04 


0.458 


0.364238E-04 


0.467 


0.434028E-04 


0.475 


0.493952E-04 


0.483 


0.539938E-04 
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0.492 


0.568843E-04 


0.508 
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